## functions
scaled_coef_cluster <- function(ModelResults, data, clusterid, subvar = FALSE, vars = NULL) {
  library(multiwayvcov)
  library(lmtest)
  library(dplyr)
  require(ggplot2)
  library(gridExtra)
  library(tidyr)
  library(stringr)
  library(purrr)
  
  cluster <- data[, clusterid]
  
  vcov_cluster = lapply(ModelResults, function(x) FUN = cluster.vcov(x, 
                                                                     cluster))
  modcoef = Map(function(x, y) coeftest(x, y), x = ModelResults, y = vcov_cluster)
  modcoef <- lapply(modcoef, function(x) FUN =round(x, digits = 4))  
  modelcoef = lapply(modcoef, function(x) FUN = x[, "Estimate"])
  modelse = lapply(modcoef, function(x) FUN = x[, "Std. Error"])
  varnames = lapply(ModelResults, function(x) FUN = names(x$coefficients))
  dfplot <- list()
  for (i in 1:length(modelcoef)){
    ad = 0.5 / modelse[[i]]
    modelcoef[[i]] = modelcoef[[i]]* ad
    modelse[[i]] = modelse[[i]]*ad 
    
    ylo95 =modelcoef[[i]] - qt(.975, nrow(data))*(modelse[[i]])
    yhi95 =modelcoef[[i]] + qt(.975, nrow(data))*(modelse[[i]])
    ylo90 =modelcoef[[i]] - qt(.95, nrow(data))*(modelse[[i]])
    yhi90 =modelcoef[[i]] + qt(.95, nrow(data))*(modelse[[i]])
    dfplot[[i]] = data.frame(varnames[[i]], modelcoef[[i]], modelse[[i]], ylo95, yhi95,
                             ylo90, yhi90,
                             Model = paste('模型',seq_along(modelcoef)[i]))
  }
  
  coefficientdata = do.call(rbind,dfplot)
  colnames(coefficientdata) <- c("Variables", "Coefficients", "S.E",
                                 "Low95CI", "High95CI", 
                                 "Low90CI", "High90CI",
                                 "Model.Name")
  df <- coefficientdata %>%
    dplyr::mutate(Variables = as.character(Variables)) %>%
    dplyr::arrange(Variables)
  
  if(subvar == TRUE){
    df <- df[grep(paste(vars, collapse = "|"), df$Variables),]
    #re-order the variables: if factor, the order of factor levels is used; if character, an alphabetical order ist used
    df$Variables <- factor(df$Variables, levels = vars)
  } 
  
  p = ggplot(df) + 
    geom_hline(yintercept = 0, lty = 2) +
    geom_linerange(aes(x = Variables, ymin = Low90CI,
                       ymax = High90CI),
                   lwd = 2, position = position_dodge(width = 1/2)) +
    geom_pointrange(aes(x = Variables, y = Coefficients, ymin = Low95CI,
                        ymax = High95CI, shape = Model.Name),
                    lwd = 1, position = position_dodge(width = 1/2),
                    fill = "WHITE") +
    coord_flip() + 
    theme_bw() + xlab("") + ylab("标准化的回归系数") + 
    theme(legend.position="bottom",
          legend.title=element_blank(),
          axis.text = element_text(size=14),
          text = element_text(size=14,family ='KaiTi', face = "bold"),
          plot.title = element_text(hjust = .5, size = 16, face = "bold"))
  
  return(p)   
} 


coef_comparison2 <- function (ModelResults, n.sim = 1000, varname, data, val1, val2, 
          clusterid) 
{
  require(arm)
  library(multiwayvcov)
  library(lmtest)
  library(gridExtra)
  library(viridis)
  library(ggridges)
  library(dplyr)
  library(ggplot2)
  library(tidyr)
  library(stringr)
  library(purrr)
  cluster <- data[, clusterid]
  vcov_cluster = lapply(ModelResults, function(x) FUN = cluster.vcov(x, 
                                                                     cluster))
  modSumm = mapply(function(x, y) coeftest(x, y), x = ModelResults, 
                   y = vcov_cluster)
  noModels = length(modSumm)
  model_names = paste("Model", 1:noModels)
  set.seed(12345)
  sim = Map(function(x, y) mvrnorm(n = n.sim, coef(x), y), 
            ModelResults, vcov_cluster)
  fd <- list()
  for (i in 1:length(ModelResults)) {
    X1 <- model.matrix(ModelResults[[i]])
    X2 <- model.matrix(ModelResults[[i]])
    X = model.matrix(ModelResults[[i]])
    value = as.numeric(quantile(X[, varname[i]], c(val1, val2)))
    X1[, varname[i]] = value[1]
    X2[, varname[i]] = value[2]
    fd[[model_names[i]]] = apply(apply(X2, 1, function(x) plogis(sim[[i]] %*% 
                                                                   x)) - apply(X1, 1, function(x) plogis(sim[[i]] %*% 
                                                                                                           x)), 1, mean)
  }
  df_plot <- map(fd, data.frame) %>% map2_df(., names(.), ~mutate(.x, 
                                                                  id = .y))
  names(df_plot)[1] <- "value"
  p = ggplot(df_plot, aes(x = value, y = id, height = ..density.., 
                          fill = id)) +
    geom_density_ridges(col = "grey70", scale = 0.8, 
                        show.legend = F) + 
    scale_fill_viridis(discrete = TRUE) + 
    geom_vline(xintercept = 0, colour = "black", lty = 2) + 
    theme_ridges(font_size = 13, grid = F, center_axis_labels = T) + 
    theme(axis.title.y = element_blank())
  return(p)
}
